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ABSTRACT 

We have simulated 2.5 xlO 3 seconds of the late evolution of a 23M Q star with full hydrodynamic 
behavior. We present the first simulations of a multiple-shell burning epoch, including the concur- 
rent evolution and interaction of an oxygen and carbon burning shell. In addition, we have evolved 
a 3D model of the oxygen burning shell to sufficiently long times (300 seconds) to begin to assess 
the adequacy of the 2D approximation. We summarize striking new results: (1) strong interactions 
occur between active carbon and oxygen burning shells, (2) hydrodynamic wave motions in noncon- 
vective regions, generated at the convective-radiative boundaries, are energetically important in both 
2D and 3D with important consequences for compositional mixing, and (3) a spectrum of mixed p- 
and g-modes are unambiguously identified with corresponding adiabatic waves in these computational 
domains. We find that 2D convective motions are exaggerated relative to 3D because of vortex insta- 
bility in 3D. We discuss the implications for supernova progenitor evolution and symmetry breaking 
in core collapse. 

Subject headings: hydrodynamics, turbulence, stars: interiors, core collapse 



1. INTRODUCTION 

Numerical simulations of stellar evolution are gener- 
ally based upon restrictive assumptions regarding dy- 
namics (e.g., hydrostatic balance and mixing-length con- 
vection), because the dynamic timescales are so much 
shorter than the nuclear burning timescale. Neutrino 
cooling accelerates the last burni ng stages so that direct 
dynamic simulation i s feasible ijBazan fc Arnettl Il998t 
lAsida fc ArnettH 200rJ), at least for the oxygen burning 
shell. We are extending this work to longer evolutionary 
times, larger computational domains, and three dimen- 
sional flow (3D). In this letter, we summarize new re- 
sults with a discussion of the hydrodynamics underlying 
important symmetry breaking and compositional mixing 
processes which may significantly affect progenitor and 
core-collapse supernova models. A detailed discussion 
of these results will appear separately; we indicate the 
wider implications here. 

2. A DOUBLE SHELL MODEL: ACTIVE OXYGEN AND 
CARBON BURNING 

Previously we have evolved a 23M Q model with the 
one-dimensional TYCHO code to a point where oxy- 
gen and carbon are burning in concentric convective 
shells whi ch overlay a silicon-rich core; details may be 
found in (jYoung fc Arnettl l2005|) . The resultant stel- 
lar structure is presented in Figure For subse- 
quent hydrodynamic evolution we use PROMPI, a ver- 
sion of the PROMETHEUS dire ct Eulerian PPM code 
(jFrvxell. Muller. fc Arnettl H.989) that has been ported 
to multi-processor computing systems via domain de- 
composition. Our multi-dimensional calculations use the 
same physics as the one-dimensional TYCHO code, in- 
cluding nuclear reaction rates, equation of state, and ra- 
diative opacities. We use a 25 nucleus reaction network 
in these models, tuned to capture the oxygen and car- 
bon burning energy generation rates to within 1% of the 
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177 species version used to evolve the ID model. The 25 
nucleus network contains electrons, neutrons, protons, 
4 He, 12 C, 16 0, 20 Ne, 23 Na, 24 Mg, 28 Si, 31 P, 32 S, 34 S, 
35 C1, 36 Ar, 38 Ar, 39 K, 40 Ca, 42 Ca, 44 Ti, 46 Ti, 48 Cr, 50 Cr, 
52 Fe, 54 Fe, 56 Ni, and all significant strong and weak in- 
teraction links. The reaction rate s , incl uding 12 C(a,7), 
are from iRauscher fc Thielemannl ((2000) . 

A two dimensional model has been calculated on a 90° 
wedge which is embedded in the equatorial plane of a 
spherical coordinate system and has radial limits which 
encompass both the oxygen and carbon burning convec- 
tive shells. Table 1 lists some additional details of the 
simulated model. A three dimensional model including 
just the oxygen shell and bounding stable layers is par- 
tially evolved at present (300 seconds). 

3. RESULTS 

Flow Topology with Two Burning Shells. Following the 
readjustment of the outer boundary due to small incon- 
sistencies in the initial ID model, a quasi-steady state 
flow develops, shown in FigEl The top half of the figure 
shows the velocity magnitude; the lower shows energy 
generation. Velocities are significant even in the non- 
convective regions, but have different morphology. The 
convective regions have round patterns (vortices) with 
occasional plumes, while the nonconvective regions have 
flattened patterns (mostly g-modes) . The flow fluctuates 
strongly. New fuel is ingested from above; the oxygen 
flame shows "feathery" features corresponding to such 
fuel-rich m atter flashing as it descends. This was prev i- 
ously seen HBazan fc Amettll 99St[A^ fc Arnettl2n0rjh 
A new feature appears in the movie version of Fig. 2, 
which shows a pronounced, low order distortion of the 
comoving coordinate, squashing and expanding the ap- 
parent circles on which carbon burning proceeds. This 
is due to the coupling of the two shells by waves in the 
nonconvective region between them. This behavior seems 
robust; we expect it to persist, so that at core collapse 
this part of the star (at least) will have significant non- 
spherical distortion. 
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Fig. 1. — The 23M Q TYCHO model used as initial conditions 
for the hydrodynamic simulation is shown here, (top) Density and 
temperature profiles are shown illustrating the large number of 
scale heights simulated as well as the complex entropy structure 
due to the burning shells, (middle) The mixing length velocities are 
shown and delineate the extents of the carbon and oxygen burning 
convection zones. The dashed vertical lines mark the boundaries 
of the simulation domain, (bottom) The onion skin compositional 
layering of the model is shown here for the four most abundant 
species. 



Fig. 2. — The magnitude of the flow velocity (top) and the net 
energy generation rate, e ne t = t nuc + e„ (bottom), is shown for a 
snapshot of the simulation which includes both an oxygen and car- 
bon burning convection zone. The carbon burning convective shell 
extends to a radius of ~4.5xl0 9 cm while the figure is truncated 
at a radius of ~2.9xl0 9 cm for clarity. A weak silicon-burning 
convection zone develops at the inner edge of the grid due to a 
small boundary-zone entropy error which accumulates during the 
course of the calculation. 



TABLE 1 
Model Parameters 
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Quantity 


Value 


Stellar mass (Mq) 


23 


Stellar age (yr) 


2.3xl0 6 


Oxygen shell convective timescale a (s) . . 


~10 2 


Carbon shell convective timescale a (s) . . 


~10 3 


Hydro simulation time (s) 


2.5xl0 3 


Inner, outer grid radius (10 9 cm) 


0.3, 5.0 


Pressure scale heights across domain . . . 


~ 9 


Angular extent of grid (rad) 


tt/2 


Grid zoning, iirXn^xnj 


800x320x1 


Numer of timesteps 


~1.5xl0 6 



a These convective timescales are based on the mixing 
length theory velocities. 



Stiffness and the Source of Density Perturbations. An- 
other asymme try, nonspherica l densi t y perturbations, 
was fo und by iBazan k, ArnettJ \ 1998(1 : lAsida fe ArnettJ 
(2000). The fluctuations in density and temperature, 
presented in the top panel of Figure as root mean 
square deviations from an angular mean, reach values 
as large as ~10% and are localized at the nonconvective 
region just beyond the convective boundary (top panel). 
The fluctuations are coincident with regions where the 
buoyancy frequency, vb, is large, which can be seen in 
the bottom panel of Figure 03 Here , v% is a me asure 
of the "stiffness" of the stratification l(Turnerlll97. c iD . and 
is proportional to the restoring buoyancy force on per- 
turbed stellar matter, 



where g is the gravity, and the term in parentheses is the 
difference between the fractional density gradient of the 
stellar structure and the fractional density change due to 
a radial (adiabatic) Lagrangian displacement. Regions 
where vb is zero are unstable to convective motions. The 
spikes in vb in our model are due to steep, stabilizing 
composition gradients which separate fuel from ash and 
lead to sharp gradients in density. Convection excites 
wave motions in the adjacent stable layers which give 
rise to the density perturbations. Similar internal wave 
phenomena can be observed in laboratory ice- water con- 
vection experiments where the largest temperature fluc- 
tuations are measured immediately above the convecting 
layer where the buoyancy frequency is large ( Townsend 
1966), highlighting the generality of this phenomenon. 

Resonant Modes. The underlying stellar structure de- 
termines the set of discrete resonant modes that can be 
excited. The narrow stable layers which bound the con- 
vective shells in our simulation, including the (truncated) 
core layer, are isolated enough from other wave propaga- 
tion regions to act as resonating cavities. These modes 
are the deeper, interior counterparts to the modes ob- 
served in helio and asteroseismology studies of milder 
evolutionary stages. 

Each mode can be uniquely identified by its horizontal 
wavenumber index, I, and its oscillation frequency, lu. We 
identify excited modes in our simulation by isolating spa- 
tial and time components of the motion through Fourier 
transforms. In Figure 0] we present a power spectrum at 
each radius in the simulation for motions with I = 4, the 
largest horizontal scale that can fit into the 90° wedge 
simulated. A direct comparison between modes identi- 
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Fig. 3. — (top) Nonspherical density and temperature per- 
turbations (root mean square fluctuation about the angular mean 
divided by the angular mean) are shown for the inner 1.2xl0 9 cm 
of the simulation domain, (bottom) In this propagation diagram 
(showing oscillation frequency versus stellar radius) the gravity 
wave and (1=4) acoustic wave propagation zones are indicated by 
the horizontal and cross-hatched shading, respectively. The grav- 
ity wave cavity is bounded above by the buoyancy frequency (solid 
line), while the acoustic cavity is bounded below by the Lamb fre- 
quency (dashed line). 

fied in the simulation and those calculated from the lin- 
earized (non-radia l) wave equation of stellar oscillations 
ljUnno et alJll989]) is presented in the right hand panel 
for two modes with significant power. Although the sim- 
ulation data has additional features, including "noise" 
in the convection zones, the mode shapes in both veloc- 
ity components are strikingly similar between simulation 
and the wave equation; the identification is unambigu- 
ous. Gravity waves evanesce (exponentially attenuate) 
beyond the boundaries of the stable layers but still con- 
tain significant power in the convection zones. Acoustic 
waves are free to propagate in the acoustic cavity which 
overlaps the carbon burning convection zone. 

During the late evolutionary epoch simulated here, the 
g-mode and p-mode propagation zones are not widely 
separated in radi us, allowing wave modes of mixed char- 
acter to couple l)Unno et alJl!989|) . The modes in the 
acoustic cavity are trapped by the boundary conditions 
of the calculation but would otherwise propagate into 
the stellar envelope where they would deposit their en- 
ergy through radiative damping, providing an additional 
channel for energy transport out of the burning region. 

The good agreement of the numerical modes with the 
analytic modes indicates that our numerical procedures 
give an excellent representation of the hydrodynamics 
of waves, even of very low mach number. We note that 
anelastic codes will not reproduce the p-mode and mixed 
mode waves properly, as we have ascertained by direct 
integration of the anelastic wave equation. 

Kinetic Energies and Wave Induced Mixing. During 
the simulations, convective motions excite waves and 
build up significant kinetic energy in nonconvective re- 
gions. The integrated kinetic energies are 2.1 x 10 45 ergs 
in the inner nonconvective region, 5.2 x 10 47 ergs in the 
oxygen burning convective shell, 9.7 x 10 45 ergs in the in- 



termediate nonconvective region, and 6.9 x 10 ergs in 
the carbon burning shell. The kinetic energy is small in 
the outer stable region, but is still increasing by the end 
of the simulation. 

In our simulations, the importance of the excited g- 
modes in the stable layers lies primarily in the role they 
play in mediating mixing at the convective boundaries 
and in the stably stratified layers. We identify a cycle in 
which kinetic energy builds up in the stable layer until 
the underlying wave modes reach non-linear amplitudes, 
breakdown and drive mixing. This process is analogous 
to the physical picture which underlies semiconyective 
mixing i|Stevensonlll97ft ILauger et a,ll1983t ISnruirlH 992ft 
but is driven on a hydrodynamic rather than a ther- 
mal timescale. The growth time for the fastest growing 
modes in our simulation is ~200 seconds and leads to 
an average migration speed of outer oxygen shell bound- 
ary of ~4xl0 4 cm s — 1 , entraining mass into the con- 
vection zone at a rate of ^1CP 4 M Q s~ 4 , significantly 
affecting the evolution. Identifying the spectrum of ex- 
cited modes in numerical simulations, including ampli- 
tudes and waveforms, provides guidance for developing 
and testing a quantitative model of this mixing mecha- 
nism. An important parameter controlling the boundary 
entrainment rate is the gradient Richardson number such 
that steeper density (and composition gradients) will 
lead to lower mixing rates s o that sharp gradients are ex- 
pecte d to form and persist l)Peltierll2003l lAlexakis et alJ 
120041) . In addition to the mixing associated with wave- 
breaking, enhanced compositional diffusion can be driven 
by the presence o f the oscillatory flow setup by g-mo des 
l)Press k Rvbickilll981trKnobloch & Merrvneldlll992j) . It 
has been demonstrated that the structure of presuper- 
nova iron cores is very sensitive to how mixing is han- 
dled at convective boundaries with significant implica- 
tions for both th e explosion mechanism an d nucleosyn- 
thetic yields (e.g. iWooslev fc Weaver! [19 8 8th If the am- 
plitudes of the wave motions identified in our simulations 
remain robust to the numerical limitations (e.g., resolu- 
tion, domain size) then neglecting the mixing processes 
associated with these waves constitutes a large source of 
error in progenitor models. 

4. DISCUSSION 

Differences in 2D and 3D. While these 2D simulations 
allow us to see the interaction of carbon and oxygen 
shells, and show wave generation at convective bound- 
aries, they impose an unphysical symmetry on the prob- 
lem. Our 3D simulations have only been carried to 300 
seconds of stellar time, but show that the wave genera- 
tion is a robust result. The flow in the convection zones, 
however, are qualitatively different between the 2D and 
3D models: the 2D flow is dominated by vortices which 
span the convection zone, while the 3D flow is charac- 
terized by smaller scale plumes. Quantitatively, the am- 
plitude of the 2D convective motions are larger than the 
3D by a factor of 8, while the 3D motions are larger than 
mixing length values by a factor of 1.5. Two-dimensional 
simulations of gravitational collapse will be misleading at 
least to the extent that convective motions are important. 

Presupernova Models. Perhaps the most important 
impact that internal waves have on stellar structure in 
the late stages of massive star evolution is the degree to 
which they drive compositional mixing. In our Simula- 
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Fig. 4. — (left) The frequency power spectrum of the 1=4 component of the radial velocity (v x ) is shown as a function of radius. The 
dashed horizontal lines indicate the frequencies of two waves with significant power that are compared to linear theory in the adjacent 
panel, (right) Two 1=4 wave forms are shown. The wave forms extracted from the simulation data (dotted line) and calculated with linear 
theory (solid line) are shown for comparison. For each mode the horizontal (v y ) and radial velocity (v x ) components are presented in units 
of cm s _1 . The shaded areas correspond to stably stratified regions. Both of these modes have some p- and g-mode character. The mode 
in the bottom panel is predominately of g-mode character, while the mode in the top panel is more clearly of mixed type. 

tions internal wave modes grow to non-linear amplitudes 
and mix material at convective boundaries on a hydro- 
dynamical timescale. Given the strong dependence of 
presupernova structure on the rate at which mixing oc- 
curs at convective boundaries we see the incorporation 
of internal wave physics into stellar evolution codes as a 
neccesary refinement. 

Symmetry Breaking. Spherical symmetry in presuper- 
nova models is broken by (1) the density perturbations 
induced by turbulence within the convection zone, (2) 
the wave interactions between burning shells, and (3) 
rotationally induced distortions. The perturbations by 
waves which are trapped between the oxygen and carbon 
burning shells are correlated on large angular scales, as 
is rotation, while the turbulent perturbations have both 
a smaller scale and amplitude. Our restricted simulation 
domain filters out wave modes with I < 4, so it is likely 
that even larger scale perturbations exist in real stars. 
Symmetry breaking will seed i nstabilities in an ou tward 
propagating supernova shock ijKuranz et al J 12005). and 
in the collapsing core. The converging case has been in- 
tense ly studied for inertial confinement fusion ijLindl. .T.I 
1998). The diverging case has implicatio ns for the prob- 
lem of 56 Ni and 56 Co decay in SN1987A (Her ant fc Ben3 
IT99"TllKifonidis et a.Ul2fK^ . 



The conclusion that internal wave modes do not grow 
to large amplitudes during core collapse through nuclear 
driven overstability (Murph y et alJl2004D . is based upon 
an analysis that ignores (1) the dynamics of convective 
motion and (2) the shell-shell interactions, both of which 
are expected to become more violent as collapse is ap- 
proached. Asymmetries in core collapse have implica- 
tions for pulsar birth kicks, explos ion mechanisms, and 
for gr avitational wave generation. iTalon fc Charbonnell 
( 2005) have shown that internal gravity waves can trans- 
port angular momentum at a rate sufficient to be impor- 
tant in the evolution of solar mass stars; we suggest that 
they are important for evolution of more massive stars to 
core collapse, and for plausible prediction of the angular 
momentum distribution in that collapse. 
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